Identification of pyroptosis related subtypes and tumor microenvironment infiltration characteristics in breast cancer

Understanding the association of pyroptosis with tumor progression, prognosis and effect on immunotherapeutic response in breast cancer (BC) is limited. This study analysed forty pyroptosis-related genes to construct the pyroptosis score. Association of the pyroptosis score with the overall survival, clinical features, tumor mutation load, immune cell infiltration, and treatment sensitivity of patients with BC was analysed. Out of 983 BC samples, 304 (30.93%) had genetic alterations with the highest TP53 frequency. We identified three separate subtypes associated with pyroptosis action. These subtypes correlate with the clinicopathological characteristics, TME immune cell infiltration, and disease prognosis. Based on the expression levels of the pyroptosis genes, we divided the pyroptosis score into a high group and a low group. The immune-activated pyroptosis subtype had a higher score with a better prognosis. We also observed that the pyroptosis score correlates with the tumor mutation burden. The pyroptosis score and disease prognosis were directly proportional. A higher pyroptosis score indicated a better prognosis. Results suggest that the pyroptosis-related gene prognosis model is closely related to the immune cell infiltration of BC. The three pyroptosis subtypes associated with BC assist in accurately identifying the tumor subtype, the prognosis of immunotherapy drugs and the patient’s therapeutic response.

Survival analysis. Univariate Cox regression analysis was used to analyze the prognostic method and P < 0.05 was regarded as statistically significant. Correlation of pyroptosis genes to BC ( Table 2). The survival curve was plotted by the Kaplan-Meier.
Consistency cluster analysis and differential gene analysis. Consensus clustering software package (https:// bioco nduct or. org/) was used for clustering and typing on the basis of the co-expression of 40 pyroptosis regulatory genes. We used the "sva" package to merge GSE88770, GSE41119 and GSE42568 GEO datasets as a validation dataset, then analyzed and validated our pyroptosis subtypes. A single sample gene set enrichment analysis (ssGSEA) algorithm helped in measuring the immune cells (28 kinds of immune cells) content following a comparison of the infiltration of immune cells among different types of pyroptosis 34 . Differentially expressed www.nature.com/scientificreports/ genes (DEGs) were obtained from the TCGA dataset of the GEO dataset. Subsequently, univariate Cox analysis was performed, P < 0.05 was filtered to obtain DEGs related to patients' prognosis, and consistent cluster classification was carried out according to the prognostic DEGs to the pyroptosis subtype based on DEGs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes 35 (KEGG) were used for the analysis of DEGs.
Scoring construction and analysis of immune-related indexes of pyroptosis. Principal component analysis 36 (PCA) was used for the evaluation of the prognostic DEGs to attain a P-score (Pyroptosisscore = PCA1 + PCA2). The relationship between P-score and immune cells, TMB, different clinical features (Gender, Age, Stage, T, N, M), and PD-L1 was analyzed.
Construction of pyroptosis related prognosis model. We performed a single factor Cox regression analysis on DEG. All BC patients were randomly classified in a ratio of 1:1 into two separate groups called the training group (n = 708) and a test group (n = 708). Furthermore, we used the Lasso Cox regression algorithm, the "glmnet" R package helped in minimizing the risk of overfitting, and the risk prediction model was established by 10× cross-validation 37,38 . Multivariate Cox analysis helped in the selection of candidate genes (Supplemental Table S1) and also in the establishment of a risk score (risk score) = Σ (EXPI × coefi)), coefi while the EXPI represents the risk coefficient and expression of each respective gene. The 708 patients in the training www.nature.com/scientificreports/ group were split into two groups based on their median risk score: a low-risk group and a high-risk group. We then ran the Kaplan-Meier survival analysis and derived a receiver operating characteristic curve. Evaluate the immune status between high and low-risk groups. CIBERSORT 42 was used to measure the quantity of 22 invasive immune cells in heterogeneous samples of the low and high-risk groups in order to estimate the ratio of TICs in TME. The correlation between the P-score and 22 types of infiltrating immune cells, as well as the association between TMB and P-score was investigated in this study.

Comparison of pyroptosis-related gene
Analysis of gene mutation and drug sensitivity. The "maftools" R package was used to build the mutation annotation format (MAF) in the TCGA database for observing somatic mutations in BC patients between the high and low-risk categories 43 . In both groups, the TMB score was also measured for each BC patient. To observe the difference in the efficacy of chemotherapeutic drugs between the two groups, the "pRRophetic" package was used to calculate the half maximal inhibitory concentration (IC50) of these drugs that are widely used for treating BC.
Establishment and verification of nomograph scoring system. Using the "rms" software tool, the clinical traits and risk scores were combined to build a prediction nomogram 44 . A score was assigned to each variable in the nomograph scoring system, and then added the scores of all variables for each sample to get the total score. For 1, 3, and 5-year survival, time-dependent ROC curves were used to evaluate Nomograms. The values between the anticipated 1, 3, and 5-year survival events and the observed results were described using the nomogram calibration plot.

Ethical approval. As this work is a bioinformatics analysis, ethical approval is not required. All methods
were performed in accordance with the relevant guidelines and regulations.

Statistical analysis.
We normalized all RNA-Seq data by the ComBat function in the sva software package.
Wilcoxon rank sum test was performed to check the difference of gene expression between normal tissues and tumor tissues. The survival curve was drawn with the help of the Kaplan-Meier method, clustering classification was carried out by consensus clustering software package, and the ssGSEA algorithm helped in the evaluation of tumor-infiltrating immune cells. All statistics were completed using the R language software package (https:// www.r-proje ct). We considered P-value < 0.05 as significant.

Consent for publication.
All author knows the situation and agrees to publish.

Results
Variation and prognosis of pyroptosis regulatory genes in BC. Initially, we studied the mutation frequency of CNV, insertion, and deletion of copy number of 40 pyroptosis regulatory genes found in BC ( Fig. 1A,B). At the same time, researchers looked at the expression of 40 pyroptosis genes in malignancies and normal tissues. CASP5, CHMP4A, CHMP7, GSDMA, HMGB1, IL1A, NAIP, NLRC4 and TP53 expressions were determined to be the same (Fig. 1C). Changes in pyroptosis regulating genes were detected in 304 of 983 BC patients, with a frequency of 30.93%. Missense mutations, splice-site mutations, and nonsense mutations are the most common types of mutations. The most frequently mutated gene was TP53, which was followed by CASP8 and DHX9 (Fig. 1D). To learn more about how the pyroptosis regulatory genes interact, we developed a network diagram of survival and interaction between pyroptosis regulatory genes ( Fig. 1E). Survival analysis revealed that 29 PRGs were closely associated with prognosis. Patients with high expression of AIM2, CASP1, CASP4, CASP8, CHMP2A, CHMP4B, CHMP6, CHMP7, ELANE, GSDMD, GZMA, GZMB, IL1A, IL1B, IL18, IRF1, IRF2, NAIP, NLRC1, TP63 and ZBP1 had a better prognosis, while patients with low expression of these genes had a better prognosis. Patients with low expression of BAK1, CHMP2B, CHMP4C, CYCS, DHX9, GSDMB, GSDMC and NLRC4 had a better prognosis ( Supplementary Fig. S1). These findings reveal that the expression of pyroptosis regulating genes differ significantly between normal and malignant tissues. Simultaneously, it has been established that pyroptosis regulating genes influence the prognosis of BC patients.
Pyroptosis subtypes and gene subtypes based on regulatory genes. The TCGA-BC data set and GSE20685 data set were clustered on the basis of 40 pyroptosis regulatory genes. It is categorized into three pyroptosis subtypes based on the cumulative distribution function (CDF) curve and area under the curve (AUC) of consensus score ( Fig. 2A). The three pyroptosis regulating gene cluster groupings demonstrated statistical differences in survival (Fig. 2B). The major components of pyroptosis regulating genes were evaluated using the PCA. Three pyroptosis subtypes were discovered to be easily distinguishable (Fig. 2C). Similar results can be obtained for our validation data ( Supplementary Fig. S2A-D). We next used ssGSEA to look at the number of immune cells infiltration in BC tumor samples and examined the difference between the three pyroptosis www.nature.com/scientificreports/ subtypes. It was observed that all three pyroptosis subtypes had a considerable number of immune cells and a high level of infiltration (Fig. 2D). The link between the three pyroptosis subtypes, pyroptosis regulatory genes, and clinicopathological parameters was then investigated. The survival rate of stage I patients increased dramatically in the pyroptosis cluster C age ≤ 50 years, and the expression of the pyroptosis gene increased substantially (Fig. 2E). Pyroptosis cluster C demonstrated enrichment pathways related to immune activation, including T cell receptor signal pathway, B cell receptor signal pathway, nod like receptor signal pathway, toll-like receptor signal pathway, chemokine signal pathway, cytokine receptor interaction, and JAK/STAT signal pathway according to GSVA enrichment analysis based on the three subtypes of pyroptosis regulatory genes (Fig. 2G,H). Pyroptosis cluster B was greatly linked with immunosuppression (Fig. 2F). We identified pyroptosis cluster C as an immunoinflammatory phenotype defined by adaptive immune cell infiltration and immunological activation based on the results of the aforementioned analysis. Pyroptosis cluster B is classified as an immunosuppressive phenotype known for its immunosuppression. Pyroptosis cluster A serves as a link between Pyroptosis cluster C and B. After that, gene analysis was carried out within the pyroptosis cluster groups, 2256 DEGs obtained by univariate Cox analysis (Fig. 3A). The consistency clustering was repeated, and three gene clustering types were identified as gene cluster A, B, and C (Fig. 3B). The results of survival analysis revealed a significant variation in prognosis (P = 0.002, Fig. 3C). At the same time, researchers analyzed the expression of 40 pyroptosis genes in tumors and normal tissues. The expression of the remaining 35 genes was different, except for APIP, CHMP6, CYCS, ELANE and TP53, which had no difference in expression (Fig. 3D). DEG-based analysis revealed that these subgroups had distinct clinicopathological features. Gene cluster A has a favorable prognosis for T1, N0, and stage I cancers (Fig. 3E).

Enrichment function of DEGs and clinical correlation analysis. Through GO enrichment analysis
of DEGs, the top 5 biological processes containing T cell activation, regulation of cell-cell adhesion, leucocyte cell-cell adhesion, regulation of T cell activation, regulation of leucocyte cell-cell adhesion. The top five cell components containing external side of the plasma membrane, secret granular membrane, membrane microdomain, membrane raft, and specific granular. And the molecular functions which include cytokine receptor binding, immune receptor activity, cytokine activity, cytokine receptor activity, and cytokine binding (Fig. 3F).  www.nature.com/scientificreports/ The cytokine-cytokine receptor interaction pathway, chemokine signaling pathway, Epstein-Barr virus infection pathway were considerably enriched in KEGG pathway analysis (Fig. 3G).

Construction of pyroptosis score and its clinical significance. The pyroptosis regulatory gene has
been discovered to have a regulatory effect on the breast cancer prognosis, cytokines, and immune infiltration of breast cancer. These conclusions, however, are predicted by the BC results. Presently, they cannot predict the pattern of pyroptosis regulatory genes accurately in a single BC patient. As a result, pyroptosis-score (P-score = PCA1 + PCA2) was used to quantify the pattern of pyroptosis regulatory genes in individual BC patients, and predict the patient's treatment response and prognosis using PCA based on DEGs. The pyroptosis subtype was linked to the gene subtype's pyroptosis score, with individuals having a higher score having a better prognosis (P = 0.024, Fig. 4A). A positive prognosis was indicated by Cluster A and Cluster C as well as higher pyroptosis scores (Fig. 4B-D). Simultaneously, the chosen DEGs regulate T cell activation and regulation, as well as cytokines and chemokines, and are closely linked with the clinicopathological characteristics. Activated

CD4+ T cells, CD8+ T cells, B cells, dendritic cells, natural killer T cells, regulatory T cells, T follicular
Helper cell, and type 1 are the immune cells among which the pyroptosis score was favorably associated with T helper cells (Fig. 4E). TMB is closely linked to a patient's prognosis. The prognosis of low and medium TMB in BC patients is better according to this study (Fig. 4F). Patients with a high pyroptosis score, even if they have a high TMB, have a better overall prognosis (Fig. 4G). TMB, genotyping, and pyroptosis scores are all positively linked (R = 0.19, Fig. 4H). 84 of 119 BC patients with high pyroptosis scores had gene mutations with a frequency of 70.59% (Fig. 4I). 625 of 853 BC patients with low pyroptosis scores had gene mutations with a frequency of 73.27% (Fig. 4J). A high pyroptosis score is a consistent predictor of outcome in BC patients under the age of 50  Fig. S3A). On the other side, we discovered a link between pyroptosis score and age, gender, lymph node metastasis, stage, and size of the tumor (Supplementary Fig. S3B). Immunotherapy. Breast cancer immunotherapy now focuses primarily on PD-L1 and CTLA4. The expression of PD-L1, CTLA4, PDCD1, PDCD1LG2, and HAVCR2 was found to be significantly higher in the high pyroptosis score group (Fig. 5A-E). We discovered that when CTLA4 was positive or PD-L1 was positive, and both CTLA4 and PD-L1 were positive, the immunotherapy scores were higher than in patients with negative CTLA4 and PD-L1 (Fig. 5F-I).

Construction of prognosis prediction model of pyroptosis gene. A single variable lasso technique
and Cox regression analysis were used to examine survival-related pyroptosis genes. To create a BC predictive risk model, Lasso regression analysis was used to create 16 pyroptosis genes risk model (Fig. 6A,B). The survival duration of the low-risk group was much longer than that of the high-risk group, according to a Kaplan Meier analysis (Fig. 6C-E). With a median risk score of 1.96 (Fig. 6I-K), the number of deaths in the high-risk group increased dramatically (Fig. 6L-N). In order to assess the risk model's predictive value in the BC cohort. The risk score ROC curve for 1, 3, and 5-year survival time (Fig. 6F-H) was further examined, indicating that it has high sensitivity and specificity for survival prediction. Meanwhile, the risk model's expression of 16 genes was assessed (Fig. 6O-Q). www.nature.com/scientificreports/ Construct a nomogram of pyroptosis risk score. A predictive nomogram incorporating risk score and clinicopathological aspects was developed to predict the prognosis of BC patients based on the difference in risk score between distinct clinicopathological variables (Fig. 7A). The calibration curve approximates the diagonal, indicating that in our nomogram, 1, 3, and 5-year OS have a good predictive capacity in our nomogram (Fig. 7B). The risk score and nomogram are effective predictors (Fig. 7C). The low-risk group was shown to be more susceptible to cisplatin and docetaxel. The low-risk group was more sensitive to both Cisplatin and Docetaxel (Fig. 7D,E). We determined through univariate and multivariate Cox regression analysis that the risk score was an independent prognostic factor influencing breast cancer patients (Supplement Fig. S5A,B).

Correlation between risk score and characteristics of tumor immune microenvironment. It
demonstrates that the prediction model's risk score is closely linked to immunity. Further analysis revealed that the risk score was positively correlated with M0 macrophages, M2 macrophages, mast cell activation, and NK cell resting content, and negatively correlated with B cells naive, dendritic cells resting, CD4 memory T cell activation, T cells CD8, T cells regulatory (Tregs) and M1 macrophages, implying that immune cell infiltration in the high-risk group was reduced, resulting in a decline of immune function (Fig. 8A). The relationship between 16 genes in the model and the number of immune cells was evaluated. It was discovered that these 16 genes were   www.nature.com/scientificreports/ highly linked to the majority of immune cells (Fig. 8B). The analysis indicated that a low TME score is strongly associated with a high immune score, while a high TME score is closely related to a high matrix score in order to further investigate if the risk score may be utilized as an immune index (Fig. 8C). And in the high-risk group TMB was higher (Fig. 8D), TMB also positively connected with the high-risk score (Fig. 8E), implying that the high-risk group is more likely the failure immunotherapy.

Discussion
Breast cancer is the major cause of cancer-related fatalities in women (approximately 15% of all cancer-related deaths in women) 1 . In addition to surgery, targeted therapy and chemotherapy are frequently used to control/ shrink bigger tumors and lower the risk of recurrence and metastasis 45 . After treatment, most tumors trigger programmed cell death 46 , and associated cell death killing breast cancer cells. Caspase-3 activation by chemotherapeutic medicines causes secondary necrosis/pyroptosis of cancer and normal cells and plays a significant role in cancer chemotherapy 47 . Pyroptosis is the formation of pores on the plasma membrane, leading to the destruction of the cell's permeability barrier and subsequent cell enlargement 48 . Active caspase-3 cleaves GSDME to create the N-terminal fragment of GSDME (GSDME-NT) when apoptosis commences. GSDME-NT will translocate and perforate, causing pyroptosis 49 . CNV is a structural variation that accounts for 4.8% to 9.5% of human genome diversity 50 . Some CNVs in TNBC can indicate poor prognosis and can act as prognostic markers, and they may be connected to lymph node metastasis 51 . Based on the TCGA cohort and GEO dataset, we initially looked at pyroptosis regulating gene mutations in BC. BC demonstrated a unique somatic mutation due to its heterogeneity, according to the findings. The pyroptosis regulatory genes of CNV had varying degrees of deletion and amplification and the mutation rate of somatic cells was as high as 30.93%, showing that pyroptosis regulatory gene mutations play an essential role in breast cancer.
To validate the accuracy of our findings using a single dataset, we carried out further association analysis using multiple other GEO datasets by cluster analyses. These datasets with different molecular features are combined to achieve improved normalization 52 . Association outcomes reflect the differences between single and grouped/cluster dataset analysis. Clustering analysis of each subtype offers concentrated groups and amplifies the molecular typing differences. Clustering analysis also improved the accuracy of tumor subtype classification 53 .
Increasing data suggest that tumor pyroptosis is related to tumor formation and progression. Pyroptosis has been to slow the growth of lung cancer tumors 54 , gastric cancer 55 , and colorectal cancer 56 . Pyroptosis can activate the innate immune system, inhibit the development of tumor cells by changing TME, and even directly kill tumor cells 39 . DEGs screened according to three pyroptosis subtypes were shown to be involved in T cell activation and cytokine interaction in this validation, indicating that breast cancer is intimately linked to inflammation and immunological modulation. Inflammatory bodies and IL-1 have been linked in recent studies showing that they play a vital function in promoting breast cancer growth and metastasis. Tumor growth is linked to an elevated level of IL-1 in the tumor microenvironment in mice mammary tumor models and human breast cancer tissues, which increases the infiltration of myeloid cells like tumor-associated macrophages (TAMs) and myeloid-derived suppressor cells 57 .
This study discovered that there was a clear enrichment of immune cells in Cluster C, as well as enrichment associated with immune activation, implying that localized death may play a role in breast cancer immune regulation. We also developed a pyroptosis score quantitative approach to identify different pyroptosis regulatory gene subtypes and serve as a guide for individual evaluation and treatment choices in this investigation. The immune activation pyroptosis pattern had a higher score and a better prognosis, according to the findings. TMB is one of the newest biomarkers in the field of cancer immunotherapy 58 , TMB is considerably greater in ER-negative BC individuals, particularly in TNBC patients 59 . The pyroptosis score is linked to the number of TMB in this study. The better the prognosis of patients with high TMB, the higher the pyroptosis score, showing that the pyroptosis score can be utilized as an independent prognostic marker.
According to the results of the analysis, the survival prognosis of the high pyroptosis score in our model was better than the low score, it was found that in the pyroptosis model of gastric cancer 60 , melanoma 61 , head and neck squamous cell carcinoma 62 , patients with high scorch death scores also had a good prognosis, and the immunotherapy effect was better than that of the low pyroptosis score group, which was consistent with our research results. However, the prognosis of low scores in low-grade gliomas 63 indicates that the scorch death score may play an important role in different tumors, which can better predict the TME status of tumors.
It has been found that PD-L1 is significantly more expressed in breast cancer tissues, especially in triplenegative breast cancer than in normal breast tissue, and the safety and efficacy of PD-L1 inhibitor pembrolizumab in triple-negative breast cancer, hormone receptor positive, HER2-negative, local recurrence or metastatic breast cancer are significantly enhanced 64,65 .
In this study, the expression of PD-L1, CTLA4, PDCD1, PDCD1LG2 and HAVCR2 in high focus death score was found to increase, and combined with the results of the literature, it can be shown that immunotherapy in patients with PD-1/PD-L1-positive breast cancer has significant therapeutic advantages and clinical efficacy. Similar therapeutic effects were also found in lung cancer 66 and melanoma 67 , consistent with our findings. Clinical trials have shown that tumor cells with higher levels of TMB are more easily recognized by the immune system and therefore have a stronger immune response to immune checkpoint inhibitors. If the tumor mutation burden is greater, there may be a good response to immunotherapy drugs (PD-1/PD-L1 inhibitors) Nivolumab, Pembrolizumab and Atezolizumab 68,69 .
Pyroptosis enhances immune activation and function, resulting in tumor clearance. Furthermore, tumor cells can activate pyroptosis in a variety of ways, and some immune cells can directly generate it, implying that pyroptosis is implicated in the positive feedback control of anti-tumor immunity 70 . In BALB/c mice treated with NP-Gsdma3 and Phe-BF3, the number of CD4+, CD8+, natural killer (NK), and M1 macrophages increased. www.nature.com/scientificreports/ Monocytes, neutrophils, myeloid-derived suppressor cells, and M2 macrophages all decreased, implying that pyroptosis may play a role in tumor immune control 71 . GSDMB-mediated pyroptosis can function downstream of GZMA, and cytotoxic lymphocytes can transmit GZMA to GSDMB-expressing cancer cells, enhancing antitumor immunity 72 . This study was a retrospective analysis using information from the database. Selective bias could skew the results; therefore, more data from BC patients undergoing immunotherapy is needed to confirm the study's conclusions. Clinical data such as surgery, neoadjuvant chemotherapy, radiation, and chemotherapy are not studied, which could alter the immune response and pyroptosis prognosis.

Conclusion
In this investigation, we genotyped BC samples using 40 pyroptosis genes to assess the effect on tumour immune matrix milieu, clinicopathological characteristics, and prognosis in BC patients. The therapeutic benefits of different subtypes and immunotherapy were investigated using a pyroptosis prognostic model. This study adds to our knowledge of the regulatory role of pyroptosis genes in BC. Findings also provide a valuable reference for guiding personalised immunotherapy and BC prognosis.

Data availability
The datasets used and/or analyzed during the current study available from the corresponding author on reasonable request.